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Abstract. We discuss the possibility that astrophysical accretion disks are dynamically unstable to non- 
axisymmetric disturbances with characteristic scales much smaller than the vertical scale height. The instability 
is studied using three methods: one based on the energy integral, which allows the determination of a sufficient 
condition of stability, one using a WKB approach, which allows the determination of the necessary and sufficient 
condition for instability and a last one by numerical solution. This linear instability occurs in any inviscid stably 
stratified differential rotating fluid for rigid, stress-free or periodic boundary conditions, provided the angular 
velocity f2 decreases outwards with radius r. At not too small stratification, its growth rate is a fraction of Q. The 
influence of viscous dissipation and thermal diffusivity on the instability is studied numerically, with emphasis on 
the case when dlnfi/dlnr = —3/2 (Keplerian case). Strong stratification and large diffusivity are found to have 
a stabilizing effect. The corresponding critical stratification and Reynolds number for the onset of the instability 
in a typical disk are derived. We propose that the spontaneous generation of these linear modes is the source of 
turbulence in disks, especially in weakly ionized disks. 
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^ • 1. Introduction 

b " 

i The simplest model of an accretion disk is that of a barotropic, axisymmetric rotating shear flow in hydrostatic 
vertical equilibrium, with a Keplerian velocity law. Realistic disks are also subject to baroclinic effects (when the 
rotation departs from cylindrical, i.e. when O varies also with the axial coordinate z) and to a vertical stratification, 
induced either by the hydrostatic state or via the illumination of the surface of the disk due to the central object. If 
the stratification is unstable, it leads to turbulence via convective instability. When the stratification is stable, it is 
generally ignored or thought to be unimportant in the stability analysis, under the rationale that it can only stabilize 
the flow. Ignoring the stratification and baroclinic effects makes the accretion disk look like a simple differentially 
rotating shear flow, with an azimuthal keplerian angular velocity profile f2(r) oc r~ 3 / 2 . Its linear stability with respect 
to axisymmetric disturbances is governed by the Rayleigh criterion in the inviscid limit: 

d(r 2 Q) 2 

-i '— > for stability. (1) 

dr 

Flows obeying this criterion are called centrifugally stable. Keplerian flow, in which angular momentum increases 
outwards, falls into this category. Yet, there is observational evidence that astrophysical (putatively Keplerian) disks 
are turbulent (c/. Hersant et al. 2003), and thus that a source of instability exists in these flows. 

Leaving aside baroclinic effects, various mechanisms have been found able to destabilize centrifugally stable flows. 
They may or may not apply to astrophysical disks. 
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i) Centrifugally stable flows can experience a globally subcritical bifurcation (Dauchot & Manneville 1997), induced 
by finite amplitude disturbances involving non-linear mechanisms not captured by the Rayleigh criterion (Dubrulle 
1993). The transition threshold in this case is related to the amplitude of the external disturbance, as typically observed 
in a plane shear flow (Dauchot & Daviaud 1994 ). 

Taylor-Couette experiments, with fluid sheared between two concentric rotating cylinders in the centrifugally stable 
regime, have revealed such transition. When the inner cylinder is at rest, the data of Wendt (1933) and Taylor (1936), 
re-analyzed by Richard and Zahn (1999), show that for the typical amplitude of the intrinsic disturbances of the 
experimental devices, turbulence subcritically sets in as soon as R = Ud/v > , where U is the relative velocity 
between the walls, d, the radial extent of the flow (the gap) and v the viscosity. 

For small curvature, the threshold is essentially independent of the gap/radius ratio and in the limit of plane 
Couette flow, ~ 1500, a value in agreement with the minimal Reynolds number for which turbulence can be 
induced in plane Couette flow. For large curvature, R™ 1 scales with the square of the gap/radius ratio. 

Finally, novel laboratory experiments were performed recently (Richard 2001), to refine the present analysis and 
to explore regimes with co-rotating cylinders. For the first time, the hysteretic behavior of the transition with the 
inner cylinder at rest has been described, a signature of subcriticality and turbulent regimes were detected when the 
angular velocity decreases outward, in the centrifugally stable region. 

ii) Centrifugally stable flows can also be destabilized by compressibility effects via non-axisymmetric instabilities 
(Papaloizou & Pringle 1984). This instability occurs independently of boundary conditions as long as dln£l/d\nr < 
— \/2> (Papaloizou & Pringle, 1985; Glatzel, 1987). However, the existence of such an instability for Keplerian disks 
requires the presence of at least one sharp edge (Goldreich & Narayan, 1985). This condition may be unrealistic in 
standard accretion disks. 

iii) Centrifugally stable flows can be further destabilized by adding another restoring force, which acts as a catalyzer. 
A first example is a vertical magnetic field (Velikhov 1959; Chandrasekhar 1960), which renders such flows unstable 
provided they are anticyclonic, i.e. if d(fl) 2 /dr < 0. The application of this mechanism to disks was first discussed by 
Balbus and Hawley (1991); they showed that the stratification of the disk does not modify the result, and that the 
maximal growth rate of instability in that case is of the order of the angular rotation velocity in the inviscid limit. The 
instability of a rotating flow subject to vertical magnetic field includes a surprising paradox: experimentally, it was 
found that in liquid metals, in the centrifugally unstable case, the magnetic field inhibits instabilities, and thus, has a 
stabilizing influence (Chandrasekhar 1960). In the centrifugally stable case, it creates an instability. Summarizing, it 
appears that a stabilizing factor acts upon a stable flow so as to generate instability. 

iv) A similar behavior was observed in rotating flows subject to a vertical stable stratification (Whithjack & Chen 
1974, Boubnov and Hopfinger, 1995). In the centrifugally unstable case, the stratification enhances the stability of 
the flow and tends to increase the critical Reynolds number of the inner cylinder (Fig.^l. In the centrifugally stable 
regime, the experimental stability curve crosses the critical line for the Rayleigh criterion, and turbulence sets in via a 
non-axisymmetric mechanism. This instability is present in the small gap and wide gap regime, showing that curvature 
effect do not play a role in the instability mechanism. A possible theoretical explanation for this instability has been 
given recently by Molemaker et al. (2001) and Yavneh et al (2001). They performed an analytical and numerical 
stability analysis of a rotating flow in the presence of a stable vertical temperature gradient, and discovered the 
existence of a linear non-axisymmetric instability for all anticyclonically sheared flows. In their work, they interpret 
this instability as being caused by the interaction of two edge modes through a mechanism of arrest and phase locking 
along the boundaries by the mean shear flow. 

We note that in astrophysical context, the stability of plane rotating shear flow has often been tackled through the 
so-called shearing sheet transformation introduced by Goldreich and Lynden-Bell (1965). In this approximation, the 
disk structure is equivalent to a rotating plane Couette flow, with shear being —3/2 of the rotation. The advantage 
of this kind of analysis is the possibility to investigate stability, by using a special class of disturbances for which the 
stability analysis of this flow can be reduced to study of ordinary differential equations. These modes, first introduced 
by Lord Kelvin in 1887, have the shape u(x,y, z,t, k, (3, k z ) = u(t) exp(i(k — Sj3t)x + i/3y + ik z z). Here S is the 
Couette shear rate, with the flow velocity in the y-direction varying with x. Individual non-axisymmetric such modes 
have been found to exhibit algebraic transient growth under various conditions, such as compressibility (Dubrulle 
& Knobloch 1992) or vertical stable and unstable stratification (Knobloch 1984; Korycansky 1992). However, any 
amount of dissipation causes ultimate decay of these modes, due to shear winding of the azimuthal wave-numbers. 
From a theoretical point of view, it can be shown that these modes are associated with the non-normality of the 
operator governing the perturbation dynamics. They may be the basis of a noise amplification mechanism, resulting 
in significant angular momentum flux (Ioannou & Kakouris 2001). However, these modes do not form a complete 
basis (for example, they cannot describe perturbations with periodic boundary conditions in x), and their study alone 
is not sufficient to determine the stability properties of the flow, as recognized by Korycansky (1992). A complete 
stability analysis requires solving the "initial condition" problem, looking at the stability properties of an infinite 
superposition of sheared modes, moving with the flow. However, even in the case of non-stratified rotating shear flow, 
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Fig. 1. Neutral-stability curves in a stratified Taylor-Couette flow experiment for different stratifications, in the domain 
(i?2,-Ri) where R2 and R\ are the Reynolds numbers based on the gap size and the velocity at the outer or inner 
cylinder, a) In the wide gap limit (ratio of inner cylinder to outer cylinder equal to 0.2), experiment of Whithjack and 
Chen (1974): N = 1.5s -1 (circles), N = 1.25s -1 (diamonds) and N = 0.88s -1 (triangles); b) in the small gap limit 
(ratio of inner cylinder to outer cylinder equal to 0.8), experiment of Boubnov and Hopfinger (1995) N — 1.21s _1 
(circles), N — 0.89s -1 (diamonds), N = 0.54s -1 (triangles) and N — 0s -1 (crosses). The continuous line is the neutral- 
stability line for centrifugal instability computed by Snyder (1968) in an un-stratificd Taylor-Couette experiment. In 
the centrifugally unstable case, one sees that the stratification enhances stability. In the centrifugally stable case, 
stratification favors instability. 



this problem is difficult to handle analytically and calls for a numerical solution (Cambon et al. 1994). Therefore, the 
shearing-sheet/sheared mode method provides only a partial answer to the stability problem. 

In the present paper, we resort to other approaches, in which boundary conditions must be fixed a priori, but 
which allow for a complete treatment of the instability. Our approach is complementary to the usual sheared mode 
approximation. We use two classical, analytically tractable approaches to reexamine the instability of a stratified 
astrophysical disk: one based on the energy method, the second based on the WKB approximation. The energy method 
is valid for a wide class of boundary conditions, namely rigid (vanishing velocity at the domain boundary), stress 
free (vanishing velocity derivative at the domain boundary) or periodic boundary conditions in the shear direction. 
Therefore, it does not apply to individual sheared modes, which satisfy none of these requirements. This method leads 
to a sufficient condition for stability valid in the inviscid limit, and for perturbations with a characteristic scale much 
smaller than the vertical scale height (Section 2.1). The influence of curvature on this instability is discussed in Section 
2.2. In Section 3, we use a method based on the WKB approximation to derive an explicit solution of the stability 
problem and exhibit unstable modes in the parameter space covering the sufficient condition for stability. This shows 
that the sufficient condition for stability is probably necessary. This theoretical study is completed in section 4 by a 
numerical study, probing instability regimes and the influence of dissipative processes on the instability. This study 
enlarges the numerical study of Yavneh et al. (2001) towards conditions more typical of astrophysical disks, namely 
Keplerian velocity profile and finite Prandtl number. In Section 5, we discuss the importance of this mechanism for 
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turbulence generation in disks, and the similarity between this instability and the instability generated by a vertical 
magnetic field. Our conclusions follow in Section 6, where the interplay between the present shear instability and its 
baroclinic counterpart is discussed. 



2. Theoretical study 

We consider a differentially rotating compressible stratified disk. In the following, we focus on perturbations with 
typical radial scales small compared with the vertical scale height of the disk. This limit provides two simplifications: 
first, it allows us to consider only the barotropic case, in which all vertical dependence of the equilibrium quantities is 
ignored. Then it allows elimination of acoustic waves and we can work in the Boussinesq approximation (Korycansky 
1992). The dynamical equations ruling the perturbation u take then the simple form 

V • u = 0, 

D t u + u • VU + 2£le 2 xu + Vp-/ig = vAu, 

Dth + u-WH = -£-Ah, ( 2 ) 
Pr 

in a frame rotating with the constant angular velocity f2 (to be chosen later). Here D t = <9f+U- V is the total derivative, 
v is the viscosity, Pr is the Prandtl number, g = VP/ p the local effective gravity, H the basic stratification, and h is 
the stratification perturbation in the Boussinesq approximation. 

The instability we are interested in is present both in rotating Couette flow (Kushner et al. 1998), i.e. for plane 
geometry, or in Taylor- Couette flow (Molemaker et al. 2001, Yavneh et al. 2001) i.e. for circular geometry. This 
suggests that curvature effects do not play a role in the instability. To clarify the presentation, we first deal with the 
simpler, plane case in Section 2.1. In the astrophysical context, this case is equivalent to the so-called shearing sheet 
approximation of Goldreich and Lynden-Bell (1972). It is also the small gap limit of the Taylor-Couette flow, and is 
relevant to many laboratory experiments. After discussion of this simple, illustrative case, we come back to the large 
gap limit (relevant to disks) in Section 2.2, where curvature effects are included. For simplicity, we also consider in this 
Section only the inviscid limit v = 0, keeping the general case of finite viscosity for numerical exploration (Section 3). 

2.1. The plane Couette case 

In this case, we assume Cartesian geometry x, y, z. The basic flow is a combination of a pure shear flow 

U = Sxe y , (3) 

with a constant rotation Q along the z axis. 

We expand the perturbation into normal modes: (u, h) — (u, h)(x) exp[z(/3y + k z z — Lot)]. Plugging this into J3J) 
and noting the velocity perturbation u = {u,v,w) and the rescaled entropy perturbation jh = 8d z H, we obtain the 
following system: 

— iou — 2Qv = —Dp 
—iav + Zu = —i(3p 
—iow + N 2 6 = —ik z p, 
—iaO = w, 

Du + ij3v + ik z w = 0. (4) 

Here, we have introduced the ^-derivative D and the Doppler-shifted frequency a = lo — (3Sx. The equation also 
involves the frequency Z = 2fl + S and the squared Brunt- Vaisala frequency in the axial direction: 

N 2 = --g z 8 z H, (5) 

7 

which is positive in the case of stable stratification. The system of equations (@J also describes a Keplerian disk in the 
shearing sheet approximation, provided £1 is the local Keplerian rotation and Z = 0/2. 

The system (0J is a coupled differential system. To investigate the basic mechanism of instability, we use JH) to 
eliminate three components in order to get a single differential system for one component, say u. For this, we first 
eliminate the pressure by using the three last equations of |Q| to get: 



.N 2 -a 2 



-i- 



{Du + i/3v) . (6) 
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We then substitute this expression into the second equation of (0J to obtain v as a function of u. This expression 
can be used in the first equation of to obtain a closed equation for u. After simplification and rearrangements, we 
obtain finally: 

(K \ BZ [K\ 9 fLZ 
D \G Du ) + a D \G) +2kz ~G U ~ U = ' (?) 

where K = a 2 — N 2 and G = (k 2 + 3 2 )a 2 — (3 2 N 2 . It is easy to check that this equation is the zero-curvature limit of 
that derived by Yavneh et al. (2001). Note that both K and G depend on x through a. 

Eq. Q is a classical generalized eigen-value problem (corresponding to a linear differential operator) which can 
be solved once the boundary conditions have been specified. The instability occurs for values of the eigen-values such 
that u>, and therefore a has a positive imaginary part. 

To derive stability conditions, we use an integral method adapted from Chandrasekhar (1960). We first reset 

equation Q into a standard shape by the change of function u = u (|r) ~ to obtain: 
= 0, 

a 2 Da aBZ „, 2 17Z G 

where a = D(\n(K/G). We then multiply JSJ by the complex conjugate u* and integrate over the whole domain. We 
obtain an integral equation, valid for both rigid (u = at boundaries), free-slip (Du = at boundaries) or periodic 
boundary conditions: 

dx\Du\ 2 + J dxE\u\ 2 = 0. (9) 

Both |-Du| 2 , and |w| 2 are positive continuous real function over the integration domain. The function E is a complex 
function of x (through er) whose real and imaginary part E r and Ei are continuous over the integration domain. By 
continuity (intermediate value theorem), there exist two points xi and x 2 such that: 

J dx \Du\ 2 + E r {xx) J dx \ii\ 2 = h + E r (x x ) I = 0, 

Ei(x 2 ) J dx\u\ 2 = 0. (10) 

In principle, a stability condition may then be found by imposing that there exist no solutions with positive u>i for 
any value of Xi and X\. In practise, the algebraic equation is of order eight in w^, and we were not able to derive 
simple explicit stability conditions. We therefore restrict ourselves to the instability condition in the limit of weak 
non-axisymmetry, by using the fact that iOi vanishes in the limit (3 = and that the instability is non-oscillatory 
uj r = (see Molemaker et al, 2001). Then, since a = —flSx + iuji, we can set a = fia and expand E as a function of 
p. To first order in /3, we find: 

e__ nz^ sn f SN V 

k 2 ' 7V2~ k 2 d 2 -N 2 + \k 2 a 2 ~N 2 ) ' 

To obtain an analytical condition, we set pe 1 ^ = k 2 a 2 — N 2 and note that finding positive cui as a function of x 2 — x\ 
is equivalent to finding p as a function of ip. In these variables, the first equation i|10fl becomes to first order in 3: 

Ii N 2 S \ 2 ( SN 2 \ 2 ( S \ 2 2 



*2/ 4n» ' m) p + { p --m C0S V (^-3)a^o. (12) 

In the limit of vanishing stratification, the last term of the l.h.s. of eq. I|12J) is negligible. The second term of the l.h.s. 
is always positive. So a (sufficient) condition which guarantees the absence of unstable solutions is that the first term 
of the l.h.s. is positive, i.e.: 

a 1 j i\j2 

m > -,W W STABILITY. (13) 

Using the energy method, we have derived a sufficient criterion for stability for periodic, stress-free or rigid boundary 
conditions. Therefore, it does not prove that flows with 5/217 < are unstable. However, in the WKB approximation of 
Section 3 and numerical stability analysis performed in Section 4, we found that all flows with 5/217 < are unstable. 
This shows that condition l|13fl is also a necessary condition for stability for these boundary conditions. 
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In the limit of vanishing stratification, we find 5/20 > as a criterion for stability, instead of the S/2VI > —1 
criterion for centrifugal stability. The stratification thus enlarges the domain of instability. However, this criterion is 
non-axisymmetric; for /3 = 0, the only non-trivial modes are stable, with the epicyclic frequency a 2 = (2ClZ)k 2 . Note 
that our criterion is in agreement with the criterion derived by Molemaker et al (2001) and Yavneh et al (2002), which 
was S/2VL < for instability. Note also that as the stratification increases, it becomes increasingly easy to satisfy 
<|13|) . So, while a weak stratification destabilizes the flow, a large stratification re-stabilizes it. We show in Section 
2.6 that when the stratification is replaced by a vertical magnetic field, a similar phenomenon occurs: weak fields are 
destabilizing, while very large fields are stabilizing. 

2.2. Influence of curvature 

We now explore how curvature effects modify the instability. For this, we assume cylindrical coordinates r, z and 
consider the stability of a general rotating shear flow 

U = rO(r)e , (14) 
Thus, the rotation frequency is noted O, and the shear rate S = rd r d. The normal-mode decomposition is 
(u, h) = (u, h)(x)e lil ^ z -^\ (15) 
The analog of in this general case is derived in Yavneh et al 2001. It is: 

D ( — rD (ru) 

i ( kz\ ,oz i\ , x 

-D[— )+2^---j M = 0, (16) 

where Z = 20 + S, K = a 2 - N 2 , G = (r 2 k 2 z + l 2 )a 2 - l 2 N 2 , a = uj - ICl and D is the derivative with respect to r. 

Following a method proposed by Dubrulle and Graner (1994), we now perform the change of variable x = ln(r/ro), 
where ro is any characteristic radius, and the change of function u = ru. Setting k z = k z r exp2cc, we may then write 
ljlU|) as: 

k„ _\ i i ' kz\ ~ 9 nz 

_UD + Oh 2 7/.-?/.= 



D x [qD x u J + -uD I — J + 2^— u -u = 0, (17) 

similar to J7J, the equation for the plane case. The main difference here is that k z , Z, Q and S are now functions of 
x. Keeping this in mind, we then repeat the same steps as in Section 2.1 until we get the equivalent of (|10fl with the 
first order expansion of E given by: 

~ 9 flZ DZ 

N 2 + k 2 (aDs^,s-4sn) n*-s)N\ 2 
k 2 & 2 -N 2 \k 2 a 2 -N 2 J 

From (|18(l . we sec that the only effects of curvature occur at small k z . Indeed, for large k z , ak z = 0(1) , we can expand 
Ijl8(l and get to first order in k z : 

E = 2fc -V - ^ N 2 (k 2 a 2 - i) + 3fc - { m - D^v 2 J ' (19) 

which exactly matches the expression obtained in the plane Couette ljll|l . We thus conclude that a condition for 
stability at large k z is (|13|l. Note that in astrophysical thin disks k z ~ (R/H)expx, and so the curvature effects are 
likely to be small everywhere except in the innermost part of the disk. In that region, however, many other physical 
processes may play an important role (like stellar magnetic field, general relativity effects, radiative processes,..) so 
that we do not think it very important to explore this peculiar effect. 

3. WKB approximation 

3.1. Method 



In the previous section, we derived stability conditions independent of the explicit shape of the solutions. In this 
Section, we derive analytical solutions using a WKB approximation (Bender & Orszag 1987, Nayfeh 1978). This 
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allows us to exhibit explicit examples of unstable cases with 5/20 < 0. The WKB approximation requires a small 
parameter. In the following, it will be convenient to choose e = S/tt and let e — > 0. We introduce the slow horizontal 
variable X = ex. In the spirit of WKB approximations, the solutions for u and p will be split into a fast varying 
complex phase @(X)/e and slowly varying functions U(X) and P(X) according to 



(u(x),p(x)) = (U(X),P(X))exp(ie(X)/e) 



(20) 



As a result, the derivative of u and p are replaced by Du — iQxU + eU x and Dp — iOxP + ePx (with fx = dxf- 
This development must be plugged into the equations of motion For this, we rewrite them as a second order 
differential system for the velocity component u and the pressure p alone 



(2VtZ-a 2 )u = -2i/3n P + iaDp 
{N 2 - a 2 ){aDu + [3Zu) = i [a 2 k 2 z - f3 2 (N 2 ~ a 2 )] p 



(21) 
(22) 



For non vanishing azimuthal wavenumbers f3 =^ the quantity a depends on the x-coordinate. We take advantage of 
the particular scaling for the growth rate: u = /30d> leading to a = /30cr with a = Q — ex. After insertion into eq. 
J2Il-(|221), we obtain: 



0(2(2 + e) -[3 z a 2 )U = -2i0P + (3a{ieP x - & X P) 



- Fr 2 /3 2 a 2 )(eaU x + (2 + e + iaS x )U) = i(3GP, 



(23) 
(24) 



where G — a 2 Fr 2 a 2 — 1 with a 2 — k 2 + f3 2 , and the Froude number being Fr = Q/N. The slowly varying functions 
U(X), P{X) are then expanded as 



U = U + eUi 



and P = P + ePx 



(25) 



Upon substituting (|25|l into the governing Eqs. Q23 p -I|24 [l one is led to a set of successive problems for (C/o , Pq), (Ui,P\) 
etc... . 

At the leading order, the governing equations for Uq, Pq are 



n(4- (3 2 <J 2 )U = -f3(2i + <7e x )Po 
Fr 2 (3 2 a 2 )(2i- aO x )U = -f3GP 

The above system admits non-trivial solutions provided that Ox is given by 

(4k 2 Fr 2 +I3 2 - a 2 Fr 2 a 2 (3 2 ) 



e 2 



A" 



(1 - (3 2 Fr 2 a 2 ) 



(26) 
(27) 



(28) 



Once O is known the slowly varying functions Uq(X) and Pq(X) are determined by the condition that there are 
non-trivial solutions for U\(X) and P\(X). At the order e the governing equations are 



0(4- p 2 a 2 )Ui+ (3{2i + (j<d x )Pi = i (3a P ox - 2W 
0(1 - Fr 2 f3 2 a 2 )(2i - aQ x )U x + PGPi = - Fr 2 (3 2 a 2 ){aU x + U ) 

The existence condition is thus 

(ipaPox - 2W )G + i{l - Fr 2 f3 2 a 2 ){<jU x + U ){2i + d& x ) = 

The following identities 

(A-(3 2 a 2 ) TT G (2i-ae x ) 
Pa = - n p(2 i + aex) U ° = - HU0 (l-FrW) = 1 

are used to transform 1|31[) in an equation for Uq alone 



. U, 



nx 



U 



Hx 
H 



1 



2i(aQ x )x 



2 i 



a aQ x (2i + ae x ) ®x 



k 2 Fr 2 



20 1 



(1- Fr 2 [3 2 a 2 ) (4- f3 2 a 2 ) 







(29) 
(30) 



(31) 



(32) 



(33) 



Equation l|33|l is a closed differential equation providing the analytical shape of the modes. Solutions corresponding 
to ct having a positive imaginary part provide the unstable modes. This requires the a priori knowledge of 0(X). This 
expression is independent of A" in a few cases discussed below. These cases provide analytical expression of unstable 
modes, as we now show. 
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3.2. Exponential solutions 

In this Section, we investigate the exponential solution, obtained when is purely imaginary. This case is readily 
obtained in at least two cases. For the particular value of the Froude number Fr = 1/2, one may check that Q 2 X = —a 2 . 
Also, in the limit: (3 — > and Fr w 0(1), considered previously by Yavneh et al. one gets Q 2 X — —4k 2 Fr 2 . These two 
cases being formally identical, we focus on the last one, which has been intensively studied by Yavneh et al. Taking 
into account the constancy of 8x, the equation (JJSJlis readily integrated to give the result 

U = C(l ± k z Fra) expfrk z FrX) (34) 

where C is a constant of integration. Upon substituting the above expression in (|2t)J) one recovers the edge modes of 
Yavneh et al. 

U=(l± k z Fra) exp[=F2fc z Fr(l + e/2)x] (35) 

written here in the asymptotic limit (1 + e) 1 ! 2 — > 1 + e/2. From their work, we thus get that unstable solutions exist 
for all anticyclonic flows S/2SI < 0, while no unstable solution exist in this limit for S/2f2 > 0. This shows that the 
criterion for stability derived in Section 2 is also necessary in this limit. 



3.3. Oscillating solutions 

The oscillatory situation will correspond to real values leading to solutions in 1)20(1 which are oscillating in space. 
They occur for example when j3 is large while Fr — > and (3Fr is of order unity so expression Q28JI is approximated 

by 



-P 2 i 



k 2 Fr 2 q 2 
(1 — (3 2 Fr 2 a 2 ) 



(36) 



According to the prescribed limit the second term in the r.h.s. of (|3(j|) has a small denominator and a numerator of 
order unity so it is dominant and &x can be written 



e 



A 



±- 



j3k z Fra 



(1-/32^2^.2)1/2' 

and after integration one gets 

6 = ±-^(l-/3 2 iVa 2 ) 1/2 . 
pFr 

Upon substitution of expression (|37l) for Ox in l|33|) , one obtains 



9 , 

dx 



U^HQx 



± 



2i 



(2i + aO x ) (3k z Fr 
The identity 

2(1 - Fr 2 (3 2 a 2 Y/ 2 = d_ 
a 3 ~ d x 



k 2 Fr 2 



a(l - Fr 2 (3 2 a 2 y/ 2 



(l-Fr 2 ^ 2 ^ 2 ) 1 / 2 



2{1- Fr 2 f3 2 a 2 ) 1 / 2 



f3 2 Fr 2 



= 



a{\- Fr 2 (3 2 a 2 y/ 2 



(37) 



(38) 



(39) 



(40) 



is used to rearrange the terms between the brackets in Eq. I|39|l . The only term which remains non- integrated is 
proportional to (k 2 z — f3 2 )Fr 2 and can be neglected when the axial and azimuthal wavenumbers have the same order 
of magnitude and Fr — » 0. Under these simplifications 



U = C 



fl/2 



(1 -Fr 2 P 2 a 2 y/' 



■ exp 



=F2i 



/3k z Fra 2 



(1-Fr 2 f3 2 a 2 y/ 2 



(41) 



Under our approximations the exponential factor in (|41|) can be set equal to unity. Introducing the notations Q = 
1 — Fr 2 j3 2 a 2 and 7 = k z /e[3Fr : the velocity component U is expressed as 



U 



rV2 r 



Q 1 / 4 



A exp (i7<3 1/2 ) + 5 exp (~hQ 1/2 



(42) 



The unknown coefficients A and B will be determined by satisfying the boundary conditions at x = ±1. Since in the 
limit (3 > 1 with l3Fr - 0(1), we have P = Q^ 2 U one can indifferently impose U = or P = at the boundaries, 
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which is of some importance in the astrophysical context where it is usually considered as more realistic to impose 
boundary conditions on the pressure. We shall define the a priori complex quantities 

Qx = Q(x = l) = l- (5 2 Fr 2 {u-e) 2 (43) 
Q 2 = Q{x = -l) = l~[3 2 Fr 2 {u, + e) 2 (44) 

with Co = ui r + iu)i . The algebraic system for A and B resulting from the boundary conditions U = or P = has 
non-trivial solutions only if the following relation is satisfied 

sin 7 (Ql /2 -C^ /2 ) = (45) 

which implies that its real and imaginary parts vanish simultaneously. We were not able to find simple general 
expressions for unstable modes, but we can exhibit analytical solutions corresponding to the neutral modes (pi = 0) 
so that Q\ and Qi are both real quantities. Moreover, we shall focus on a particular solution characterized by a fixed 
value of the frequency u) r chosen so that either Q\ = or Q 2 = 0. fn the former case we have 

0FrQ r = pFrs ± 1. (46) 

When Qi = 0, the condition [7(1) = implies that A + B = 0, and the condition U ( — 1) = reduces to 

sin( 7 Q^ /2 ) = 0, or 7 Q^ /2 = Nir (47) 

where N is an integer and Q2 = —^£U) r (3 2 Fr 2 has to be positive for its square root to be real. Substituting the 
expression 1461) in the condition 1471) gives the relationship 

between the combination b — [3Fre and the axial wavenumber k z for a solution U to exist in the form 

U=^- A M{k z /b)Q^ 2 ] (49) 

with Q = —6(1 — x)[2 + b(l — x)] when —1 < b < 0. Such solutions have been plotted in Figs, la and lb for b = —1/5 
obtained with the value Nir/kz — 4. Thus varying the value of N allows us to change the value of k z — N(n/4) and 
as a consequence the number of zeros of U. The Figs.|21and ?? corresponds to N = 6 and N = 10 respectively, and as 
expected the number of zeros is proportional to TV. The Figs. 0] and ?? drawn for b = —1/10 obtained with the value 
Nw/k z = 6, correspond respectively to N = 6 and N = 10. It must be noticed that when the instability occurs via a 
Hopf bifurcation (uj r ^ 0) the eigenfunctions break the symmetry about x = 0, the amplitude being larger near one 
side of the interval, here x = —1. This has already been discussed by Knobloch (1996) who showed that system with 
0(2) symmetry can exhibit either a steady state bifurcation or a Hopf bifurcation. 



4. Numerical study 

4.1. Method 

The condition for stability derived in Section 2 is valid for a wide class of boundary conditions, but only in the inviscid 
limit. The influence of viscosity has been studied numerically by Molemaker et al. (2001) and by Yavneh et al. (2001), 
in the limit Pr — > 00. They found that a small viscosity essentially does not change the results, and only introduces a 
critical Reynolds number, above which the flow is unstable. For example, at il/N = 0.01, S/2£l = —2/3, the critical 
Reynolds number is of order 4800, for rigid boundary conditions. In a Keplerian disk, 5/217 = —3/4 and Pr is finite, so 
that the results of Yavneh et al. do not apply. We thus performed a numerical study of equation (@J in the plane case 
to investigate the range of parameters of astrophysical interest. The stability analysis was set into a classical eigenvalue 
problem via Fourier transform in the y and z direction, and discretization of the equation over collocation points of the 
Chebyshev polynomials in the x-direction, over a domain —d/2 < x < d/2. The precision of the solution is governed 
by the number of collocation points used in the computation. We used typically 25 to 50 collocation points to get 
satisfactory precision: doubling the resolution did not alter the results. Spurious normal modes were eliminated via 
interpolation, then cross-checked, at double resolution. In the following, we present results obtained using stress-free 
boundary conditions in y and z (the case of rigid boundary conditions is discussed in Yavneh et al. (2001).) 

Dimensional analysis of equation Q shows that the stability problem is controlled by four non-dimensional pa- 
rameters: the rotation number Ro = 217/5, the Ekman number Ek = v/2Vld 2 , the Froude number Fr = 17/7V and 
the Prandtl number Pr. From these numbers, one can also build another non-dimensional number of interest, the 
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Fig. 2. Neutral mode in the WKB approximation in the limit Fr = Cl/N — » and (3 — > oo for (3FrS/il 
k z = 3tt/2. 



-1/5 and 



CD 




Fig. 3. Neutral mode in the WKB approximation in the limit Fr — > and (3 — > oo for (3FrS/Vt = —1/5 and fc z = 57r/2. 



Reynolds number, defined as: i?e = IS'ld 2 /^ = Ek 1 \Ro 1 \, Using the gap d (equivalent to the typical scale of our 
perturbation, in the disk case, see Section 2) and l/2f2 as unit of length and time, we may then write Q as: 

- ia*u - v = -Dp + Ek {dl - k 2 ) u 

-ia*v + (1 + 1 /Ro)u = -i/3p + Ek(d%- fc 2 ) v 
Fr~ 2 

-ia*w H — 9 = -ik z p + Ek (d 2 x - fc 2 ) w, 
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Fig. 4. Neutral mode in the WKB approximation in the limit Fr — > and j3 — > oo for (3FrS/Cl = —1/10 and k z = ir. 
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Fig. 5. Neutral mode in the WKB approximation in the limit Fr 
k z = 5tt/3. 



and /3 -» oo for f3FrS/Q = -1/10 and 



—i<7*6 = w + 
Du + ifiv + ik z w = 0, 



Ek 



(dl - k>) e, 



(50) 



where <7* =0-/20. 

Our numerical study was focused on the determination of maximal growth rate and wavenumber, for a given 
value of the four parameters. Due to the large number of independent parameters, it is difficult to completely explore 
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the instability phase space. Moreover, we found out that for Ro < —2, the search for maximal growth rate was 
increasingly time consuming due to the exponential narrowing of the instability branch around the maximal value. 
This phenomenon had been described in Molemaker et al (2001). We thus mainly focused our analysis on the case 
Ro > —2, with special emphasis on the case Ro = —4/3 corresponding to the Keplerian case. 



4.2. Reminder 

A few results regarding special values of the parameters can already be drawn from previous studies: in the unstratified 
case Fr — > oo, Lezius and Johnson (1971) found that the neutral stability curve is given by 

Re 2 Ro(Ro + 1) = -1706, Fr = oo. (51) 

In the inviscid case Ek — 0, Molemaker et al. (2001) found that for Fr < 1, the maximal growth rate and vertical 
wavenumber scale as: 

^ = -2i?cr 1 /Mexp(2.Ro), 

k z d = — ; = (52) 

2Fr ^/Ro(Ro+l) 

In the case Fr > 1, Ek — they found that both k z d and Wi/2f2 decay linearly with 1/Fr. 



4.3. Influence of viscosity 

To study the influence of viscosity on the instability, we set Pr — oo and vary the Ekman number, for different rotation 
number and stratifications. In Figure El we show the non-dimensional optimal growth rate, vertical and azimuthal 
wavenumber as a function of the rotation number, at Ek = 10~ 4 , for two different stratifications Fr = 3.1623 and 
Fr = 0.1. The inviscid result l|52l) is also displayed for comparison. One sees that viscosity tends to decrease optimal 
growth rate and wavenumbers. Moreover, it induces a finite critical value of the rotation number below which the 
optimal growth rate becomes negative. This critical rotation number depends on the stratification (Figure 0). It 
decreases from —1.3 at Fr = 0.02 up to —2.85 at Fr = 3.1623. We note on the figure a tendency towards leveling 
at the value Ro — —3, which may indicate a critical rotation number of —3 for Ek = 10~ 4 , below which there is no 
instability, whatever the stratification. 

In Fig EI w e show the optimal growth rate and wavenumbers as functions of stratification, at Ro = —4/3, for 
different values of Ek. For comparison, the results obtained by Molemaker et al (2001) at Ro = —3/2 and Ek = 
(for different boundary condition!) are also displayed. The growth rate tends to be proportional to N for large Froude 
numbers, but becomes a fraction of the dynamical frequency f2 for Fr > 1. One sees that the viscosity (and boundary 
conditions) has little influence on the vertical wavenumber, which scales as 1/Fr, while it introduces a sharp cut-off at 
small Froude numbers in both u/j and 8. This cut-off is pushed towards larger stratification with increasing viscosity. 
The cut-off defines a viscosity-dependent critical Froude number, below which no instability is present. Its dependence 
on the Reynolds number is approximately a power law Fr c ~ i?e -1 / 2 , as shown in Fig. EI 

The variation with stratification indicates that at small stratification, the growth rate is proportional to N, while 
at large stratification, it is proportional to f2, i.e. independent of the stratification. 



4.4. Influence of Prandtl number 

For finite Prandtl number, the stratification agent (for example the temperature) undergoes a diffusive process, which 
can have a stabilizing influence. To study this effect, we fix Ro = —4/3 and Ek = 10~ 4 , and decrease the Prandtl 
number from Pr = oo to Pr = 0.1. The resulting optimal growth rate and wavenumbers are shown in Fig. 1101 One 
sees that decreasing the Prandtl number at fixed Ekman number provides qualitatively the same effects as when 
increasing the Ekman number (previous Section): as Pr is decreased, the critical Froude number is shifted towards 
higher values. We have checked that this similarity of behavior also extends to optimal wavenumbers: no incidence on 
vertical wavenumber, varying cut-off for 0. The variation of the critical Froude number with Prandtl number is shown 
in Fig. [5] The variation of Fr c with Pr is less steep than for Re, namely Fr c ~ Pr -1 / 4 . 

For illustration we also display in Fig. 1111 a typical eigenmode solution for the radial velocity and entropy pertur- 
bation for Ro = —4/3, Ek — 10~ 4 , Fr — 0.5 and two values of the Prandtl number. One sees that these eigenfunctions 
vary smoothly over the channel width and that the change with Prandtl number is rather smooth and unimportant. 
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Fig. 6. Maximal growth rate (2a), corresponding vertical wavenumber (2b) and azimuthal wavenumber (2c) as a 
function of the rotation number 2Q/S at Ek = 10~ 4 , Pr = oo and for different stratifications. The line is the result 
of Molemaker et al. (2001), obtained at Ek = 0. 
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Fig. 7. Critical rotation number vs. Froude number (Cl/N), at Ek = 10 4 . 
4.5. Influence of azimuthal wavenumber 



Up to now, we have focused only on optimal modes (those with highest growthrate) . These modes are characterized by 
a low azimuthal wavenumber. In the application to astrophysical disks, it may be interesting to look for modes with 
larger azimuthal wavenumber, which may be more realistic for very thin disks (for which (3d <~ r/H ^> 1). Finding 
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Fig. 8. Maximal growth rate (4a), corresponding vertical wavenumber (4b) and azimuthal wavenumber (4c) as a 
function of the inverse Froude number at Ro = —4/3, Pr = oo and for different viscosity. The filled symbols correspond 
to the optimal eigen mode; the open symbols correspond to a member of the "large j3 family" , computed at Ek = 10~ 5 . 
The lines refer to the theoretical formulae proposed by Molemaker et al. (2001) for the optimal mode: u>i/2fl = (3de 2Ro 
and k z d = Ro/2Fry/l + Ro. 



such modes is very time consuming using our numerical procedure (which is optimized for the detection of optimal 
growth rate). So, we focused on the case Ro = —4/3 and only explored the case Ek — 10~ 5 , Pr — oo. In that case, we 
were able to capture one member of the large /3 family whose characteristics are shown in Fig. [S] for easier comparison 
with the optimal mode. The azimuthal and vertical wavenumber of this mode are roughly twice those of the optimal 
family. Its growth rate is about 10 times smaller. The corresponding velocity profile shows a double oscillation within 
the channel, see Fig. The mode seems to exist only at rather large stratification Fr > 1. All these features are 
reminiscent of the "large /3 mode" found in the Taylor-Couette configuration by Molemaker et al. (2001). These modes 
are characterized by a "radial" wave number k > 2, thereby producing k oscillations of the eigenmode in the radial 
direction. The wavenumber roughly follows k 2 = f3 2 k 2 . Their growth rate, in the inviscid limit, scales as |il|fc _1 / 4 e 3fi:o . 
With k = 2, this gives uii/2fL — 0.008 in rather good agreement with our finding. This strongly suggests that these 
"large (3 -modes" also exist in the plane Couette case. Since their growthrate is still a significant fraction of the rotation 
period, these modes are probably also quite relevant for astrophysical disks. 

4.6. Comparison with experiments 

It is interesting to compare our results with experimental data obtained in the same conditions (small gap limit) 
by Boubnov and Hopfmger (1995). We note that in this experiment the boundary conditions are probably rigid 
(in contrast with our numerical analysis). This comparison can then be used to probe the influence of boundary 
conditions. A difficulty arises because the experiments only provide data on the critical line, on which Ro, Fr and 
Re vary simultaneously. We thus analyze the three possible planes, in Fig. 1131 In panel a), we observe a confirmation 
of the stabilizing role played by stratification: weak stratification tends to increase the range of instability, while 
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Fig. 9. Variation of Reynolds and Prandtl numbers on the neutral stability lines at Ro = —4/3 as a function of the 
Froude number. 
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Fig. 10. Maximal growth rate as a function of the inverse Froude number at Ro = —4/3, Ek = 10~ 4 and for different 
Prandtl numbers. The line is the result at Pr — oo, drawn for comparison. 

large stratification (small Froude number) tends to restrict the domain of instability. The experimental points do not 
overlap with our numerical points, since the latter are performed at constant, and slightly sma Her Ekman number. 
In panel b), we do not observe any clear dependence of the rotation number on Reynolds number (different instability 
branches may be present). At Ro = —4/3, the Keplerian value, the critical Reynolds number seems to be between 
4000 and IO 4 , well below typical Reynolds number of e.g. circumstellar disks (Hersant et al, 2003). In panel c) one 
sees the dependence of the critical Reynolds number as a function of the stratification. At small Froude numbers, 
the experimental data seem to confirm the Fr~ 2 scaling obtained numerically. At large Froude number, we observe 
another, linear scaling, of the Reynolds number with the Froude number, Re c ~ 2000Fr. This number will be used 
below in our discussion. 

5. Application to astrophysics 

5.1. Stratification in disks 

In previous Sections, we have shown that all flows with S/2Q < are linearly unstable with respect to non-axisymmetric 
perturbations in the presence of vertical stable stratification. 

To decide whether an astrophysical disk is stably stratified, on has to solve for the vertical structure of the 
disk, which itself depends on how turbulent kinetic energy is transformed into heat through viscous friction. Most 
prescriptions for this heat release lead to stable stratifications. 
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Fig. 11. Real (continuous line) and imaginary (dotted line) part of the optimal eigenfunction for the radial velocity u 
(Fig 7a) and for the scaled entropy fluctuation 9 (Fig 7b), as a function of the x-coordinate, at Ro = —4/3, Fr = 0.5, 
Ek = 1CT 4 , and for Pr = 0.1 (thin line) and 10 (thick line). 



One exception may be the central part of disks around black holes, where radiation pressure exceeds gas pressure, 
and where for this reason entropy decreases with distance from the equatorial plane, triggering thermal convection 
(Bisnovatyi-Kogan & Blinnikov 1977). But this does not take into account the illumination of the disk by hard X-rays, 
which is observed. The presence of thermal convection has also been suggested in disks around young stellar objects 
(Lin & Papaloizou 1980). However, here again a natural process leading to stable vertical stratification is the disk 
illumination by the central object (Chiang & Goldrcich 1997). 

A rough estimate of the Brunt- Vaisala frequency in stably stratified disks may be obtained in the thin disk 
approximation: 

N2 = ih [Vad _ v] = q2 (10 2 [Vad _ v] (53) 

where the actual temperature gradient d In T/d In P is determined by the opacity law; the adiabatic gradient V a d = 
(91nT/91nP) a( j takes the value 0.4 for perfect atomic or ionized gas. Typically [V a( j — V] « 0.1, and thus N/Q <~ 0.3. 
This value is confirmed by detailed calculations. 

Radiative numerical simulations of protostellar disk (D'Alessio et al. 1998) reveal that at r = 0.5 a.u. the tempera- 
ture is of the order of lOOOif at the disk scale height and 300K in the midplane, leading to N/Cl ~ 0.3. In this regime, 
the critical Reynolds number is of the order of 1000, the growthrate is about one percent of the rotation frequency and 
the typical vertical and azimuthal wavenumbers are about 0.1/d. The largest allowed d in this case is 0.1H leading to 
vertical scale and azimuthal wavenumber of the order of unity. The typical instability should then take the form of a 
spiral mode. 

We can thus conclude that, in astrophysical disks, stable stratification is the rule rather than the exception. If this 
were not the case, one would invoke thermal convection as the obvious cause for the turbulence which ensures the 
accretion on the central object. 
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Fig. 12. Real (continuous line) and imaginary (clotted line) part of a "large /?" eigenfunction for the radial velocity u, 
as a function of the ^-coordinate, at Ro = -4/3 Fr = 200 Ek = 10~ 5 , (3d = 3.3, k z d = 41, Ui/20, = 0.0013 and for 
infinite Prandtl number. 



5.2. Magnetic field vs. stratification 

Our findings are reminiscent of what occurs in the presence of an axial magnetic field. It is therefore interesting to 
examine the degree of similarity between the two instability mechanisms using the same tools. For this, we recall 
the results obtained by Chandrasekhar (1960). First, unlike the hydrodynamical instability in stratified disks, the 
magneto-rotational instability has axisymmetric unstable modes. For these modes, the sufficient stability condition is: 

S/2n >-IJ±-¥a, STABILITY (54) 
4 fcjio \V 

where FIa = y/k%nB%/4Tvp is the Alfven frequency. Comparing this condition with 1)13(1 , we see strong similarities with 
the stratified case, the Brunt- Vaissala frequency being just replaced by the Alfven frequency. For example, strong 
magnetic fields tend to inhibit instability, in the same way as strong stratification, whereas in the limit of the weak 
magnetic field S/2Q, > is a condition for stability, as for stratified, rotating shear flows. Keplerian flows, having 
S/2Q, = —3/4 do not satisfy this stability criterion, and are therefore liable to both type of instabilities. 

In a recent detailed numerical study of the magneto-rotational instability, Willis and Barenghi (2002) found that 
the critical Reynolds number at infinite Prandtl number is of the order of 20, that is one order of magnitude less than 
for the strato-rotational instability discussed here. They also noted that the variation of the critical Reynolds number 
with magnetic Prandtl number is Re c ~ 100/\/P m , similar to the law found in Fig. 5 for stratification. However, in a 
typical disk, the magnetic Prandtl number is much lower than the Prandtl number. For example, in a protoplanetary 
disk, it is less than 10 -5 (Riidiger & Zhang 2001). This means that a typical critical Reynolds number for the magneto- 
rotational instability will be greater than 30000, hence larger than the critical Reynolds number for strato-rotational 
instability. 



5.3. Barotropic vs baroclinic instability? 

In the present paper, we only focused on perturbations with typical radial scale small compared to the vertical scale 
height, leading to a barotropic description. An opposite limit would be to consider perturbations with small vertical 
scale, leading to a situation where only the vertical dependence is considered so that fl varies with the axial coordinate 
z. In this situation, one may expect new instabilities to arise, due to the vertical shear. More generally, when the 
rotation departs from cylindrical, it may induce an axisymmetric baroclinic instability which has been described 
earlier in stellar interiors by Goldreich and Schubert (1967) and Fricke (1968). The application of this instability to 
disks has been recently done by Urpin (2003), who showed that the instability is linear, and that it proceeds on a 
thermal time-scale. More generally, baroclinic instabilities occur when there is an inclination between the density and 
pressure term, breaking the vorticity conservation. As discussed by Klahr and Bodcnhcimcr (2003), this effect is likely 
to occur in any realistic disk. This means that in general, astrophysical disks are subject to both barotropic (as shown 
in this paper) and baroclinic instabilities. The interplay between these two instabilities is a fascinating subject, and 
has been widely studied in the oceanographic context. For upwelling flows it was found for example that barotropic 
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Fig. 13. Neutral-stability curves (w* = 0) in the stratified Taylor-Couette flow experiment of Boubnov and Hopfingcr 
(1995) in the small gap limit, in parameter space; a): in the plane Ro,Fr; b) in the plane Ro,Re; c) in the plane 
Re, Fr. The filled circles are the experimental data. The open circles are the numerical data obtained in the present 
paper. 

instabilities are important in the early stages of the dynamics, while baroclinic modes dominate the late stages. The 
growthrate and wavenumbers are in between the growthrate and wavenumber of each instability (Tadepalli & Fertziger 
1998). Clearly, a similar study would be welcome in the astrophysical context. 

6. Conclusion 

In this paper, we have used analytical and numerical methods to study the instability of Kcplcrian-likc flow in the 
presence of a stable vertical stratification, for perturbations obeying rigid, stress- free or periodic boundary conditions in 
the radial coordinates that lie less than one scale height apart so as to justify an incomprcssibility approximation. This 
instability is non-axisymmetric. Its growth-rate is a fraction of the rotation rate at small stratification. This 'strato- 
rotational' instability is purely hydrodynamical and does not require the presence of any magnetic field, whatever 
small. Therefore, it is especially relevant in the context of weakly ionized disks, such as the primitive solar nebula. It is 
of course also relevant in other Keplerian disks, since the critical Reynolds number to trigger this instability is of the 
order of 10 3 , which is less than the critical Reynolds number for both the magneto-rotational instability and the finite- 
amplitude hydrodynamic instability. It may also interplay with baroclinic instabilities, resulting in mixed instabilities 
similar to those observed in oceans. However the question of which instability occurs first has little relevance, partly 
because we do not know the initial conditions, but mainly because there is no reason why the instability which has the 
fastest linear growth will dominate in the fully nonlinear state. What one should rather ask is what kind of turbulence 
exists at the high Reynolds numbers which characterize astrophysical disks. A partial answer to this question may 
already be found from Taylor-Couette laboratory measurements (Dubrulle et al. 2004), which provide information 
about the hydrodynamical regime of rotating shear flow, as well as the influence of magnetic field or stratification 
on the transport properties. When applied to circumstcllar disks, this provides parameter free-predictions about 
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accretion rates and fluctuations which are in good agreement with observations (Hersant et al. 2003). For the influence 
of processes more specific to astrophysics (radiative transport, etc.) one will probably have to wait until the numerical 
simulations are able to reach Reynolds numbers of order 10 6 or more. 
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